* Settings
version 16
do "$SSDIMed/scripts/_auxiliary/_project_settings.do"


* -------------------------------------------------------------------------------------------------
* Prepare data for figures
* -------------------------------------------------------------------------------------------------

* DI entrant outcomes, by x, where x = 
*   i.covstart_year          = calendar year of Medicare entry
*   c.unemprateatapp         = unemployment rate at time of application
*   i.unemp_q20              = ventiles of unemprateatapp
*   i.age_year_covstart_fill = age in years at Medicare entry

* Sample
*   main: individuals who enter Medicare at ages 20-62, and annual Medicare spending after first (usually partial) year of Medicare coverage
*   2262: main, but individuals who enter Medicare at ages 22-62
local sample 2262

* Sample definitions:
* - main (all ages, all years)
local sample_main 1
* - ages 22-62
local sample_2262 inrange(age_year_covstart_fill, 22, 62)
* - ages 51-52
local sample_5152 inrange(age_year_covstart_fill, 51, 52)
* - first full year
local sample_ffy (rfrnc_yr == init_yr + 1)

* Which controls were included in each regression
* y = incidence
*   local controls_01 none
*   local controls_02 i.fipscounty_firstnm_g
* y = tot_pmt, died
*   local controls_01 i.years_since_covstart
*   local controls_02 i.years_since_covstart##i.fipscounty_firstnm_g
*   local controls_04 i.years_since_covstart##i.fipscounty_firstnm_g i.male_init##i.age

* DI entrant cohort size, average Medicare spending, and mortality
local clust county_mofd
foreach y in incidence_pop_age_atapp tot_pmt died_adj {
foreach x in unemp_rate_q20_atapp {
  * local y incidence_pop_age_atapp_gap1yr
  * local x unemp_rate_q20_atapp
  
  local y_short `y'
  if strpos("`y'", "incidence") > 0 local y_short incidence
  
  * Control variables
  local ctrl 02
  
  * QC: Confirm sample
  assert "`sample'" == "2262"
  
  * Load estimates
  estimates use "$SSDIMed/results/estimates/x-i.`x'/x-i.`x'_y-`y'_controls-`ctrl'_cluster-`clust'_sample-`sample'.ster"
  
  * QC: confirm controls
  noisily di "`y': `e(absvars)'"
  if "`y'" == "incidence_pop_age_atapp" assert "`e(absvars)'" == "i.fipscounty_firstnm_g"
  else                                  assert "`e(absvars)'" == "i.years_since_covstart##i.fipscounty_firstnm_g"
  
  * Calculate y effect, CI for each level of x variable (_cons plus x-var effect)
  local addvar 0, 0, 0
  local names : colfullnames e(b)
  foreach name of local names {
    if "`name'" == "_cons" continue
    *local val = subinstr("`name'", "b.`=subinstr("`x'", "i.", "", .)'", "", .)
    *local val = subinstr("`val'",  ".`=subinstr("`x'", "i.", "", .)'", "", .)
    lincom _cons + `name'
    local addvar `addvar', `"`name' + _cons"', `r(estimate)', `r(se)'
  }
  regsave, addvar(`addvar') detail(all)
  
  * which variables to keep
  drop if var == "0"
  keep if strpos(var, "+ _cons") > 0
  
  * confidence intervals for added variables
  gen double t_ci = abs(invttail(`e(df_r)', (1 - 0.95) / 2))
  qui gen double ci_lower = coef - t_ci*stderr
  qui gen double ci_upper = coef + t_ci*stderr
  rename (coef ci_lower ci_upper) (`y_short' `y_short'_ci_lo `y_short'_ci_hi)
  keep var `y_short' `y_short'_ci_lo `y_short'_ci_hi clustvar
  
  * numeric level for the x-variable (a factor)
  gen `x' = substr(var, 1, strpos(var, "."))
  replace `x' = regexr(`x', "(b\.)|(\.)", "")
  destring `x', replace
  order `x', after(var)
  
  * Units
  * tot_pmt in $1,000
  * died in deaths per 10,000
  foreach var in `y_short' `y_short'_ci_lo `y_short'_ci_hi {
    if "`y'" == "tot_pmt"        replace `var' = `var'/1000
    if strpos("`y'", "died") > 0 replace `var' = `var'*10000
  }
  
  * Save
  tempfile `x'_`y_short'
  save ``x'_`y_short'', replace
  list
}
}

* Average unemployment rate by ventile of unemployment rate
*   Average based on 1 observation per entrant
use "$SSDIMed/data/analysis/county-month-age_entry_sample-main.dta", clear
gstats tab unemp_rate_county_atapp [aw = count], by(unemp_rate_q20_atapp)
gcollapse (mean) unemp_rate_county_atapp [aw = count], by(unemp_rate_q20_atapp)
keep if !missing(unemp_rate_q20_atapp)
tempfile x_bin_means
save `x_bin_means'


* ---------------------------------------------------------------------------------------
* Graph settings
* ---------------------------------------------------------------------------------------

* Set graph style to project default settings
*	Pass as arguments: height width (default 3.5in 6.5in)
do "$SSDIMed/scripts/_auxiliary/_project_grstyle.do" 3in 3in

local red "179 0 0"
local green "75 115 47"
local blue "48 84 150"
local bluegray "126 153 180"
local brown "140 045 004"
local orange "237 125 49"
local tan "210 180 140"

* Entry
local p 2
local color`p' `brown'
grstyle set color "`color`p''": p`p'
grstyle set symbol i: p`p'
grstyle set symbolsize 2.2pt: p`p'
grstyle set lpattern solid: p`p'
grstyle set linewidth 1.8pt: p`p'
* shaded region for confidence intervals
grstyle set color "`color`p''%20": p`p'area
grstyle set color "`color`p''%20": p`p'arealine
grstyle set linewidth vvthin: p`p'area

* Medical spending
local p 3
local color`p' black
grstyle set color "`color`p''": p`p'
grstyle set symbol i: p`p'
grstyle set symbolsize 2.2pt: p`p'
grstyle set lpattern solid: p`p'
grstyle set linewidth 1.8pt: p`p'
* shaded region for confidence intervals
grstyle set color "`color`p''%20": p`p'area
grstyle set color "`color`p''%20": p`p'arealine
grstyle set linewidth vvthin: p`p'area

* Mortality
local p 4
local color`p' `green'
grstyle set color "`color`p''": p`p'
grstyle set symbol i: p`p'
grstyle set symbolsize 2.2pt: p`p'
grstyle set lpattern solid: p`p'
grstyle set linewidth 1.8pt: p`p'
* shaded region for confidence intervals
grstyle set color "`color`p''%20": p`p'area
grstyle set color "`color`p''%20": p`p'arealine
grstyle set linewidth vvthin: p`p'area

* Arrows and lines
local p 5
grstyle set color gs8: p`p'
grstyle set lpattern solid: p`p'
grstyle set linewidth vvthin: p`p'


* -------------------------------------------------------------------------------------------------
* Figure: DI entrant characteristics, by quantile of the unemployment rate
* -------------------------------------------------------------------------------------------------

* X-related variables
local x unemp_rate_q20_atapp
local x_bin_mean unemp_rate_county_atapp

* Load data
use ``x'_tot_pmt', clear
merge 1:1 var using ``x'_died_adj', assert(match) nogen noreport
merge 1:1 `x' using ``x'_incidence', assert(match using) nogen noreport
merge 1:1 `x' using `x_bin_means', assert(match) nogen noreport
list
assert inrange(`x', 1, 20)

* ------------------------------
* Panel (a): DI entrant cohort size
if 1 {

* How many major ticks?
* How much should range extend beyond ticks (as fraction of tick step size)
local y_ticks 7
local y_margin_lo 0.4
local y_margin_hi 0.4

* y1 tick lo, step size, format
local y1_tick_lo 250
local y1_step 25
local y1_fmt "%9.2g"
local y1_unit ""

* Construct y labels and scales from parameters above
forvalues y = 1/1 {
  local y`y'_tick_hi  = `y`y'_tick_lo' + `y`y'_step' * (`y_ticks' - 1)
  local y`y'_range_lo = `y`y'_tick_lo' - `y`y'_step' * `y_margin_lo'
  local y`y'_range_hi = `y`y'_tick_hi' + `y`y'_step' * `y_margin_hi'
  di "y`y'_tick_hi:  `y`y'_tick_hi'"
  di "y`y'_range_lo: `y`y'_range_lo'"
  di "y`y'_range_hi: `y`y'_range_hi'"
  
  local y`y'_label ylabel(`y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi', axis(`y') format("`y`y'_fmt'") noticks nolabels labgap(1pt)) 
  local y`y'_scale yscale(`yalt' range(`y`y'_range_lo' `y`y'_range_hi') axis(`y') noline) 
  di `"y`y'_label:  `y`y'_label'"'
  di `"y`y'_scale:  `y`y'_scale'"'
}

* x-axis settings
local xlo 1.3
local xhi 13.8
local xscale xscale(range(`xlo' `xhi'))
local x_axis xlabel(3(1)13) `xscale' xtitle("County unemployment at application", margin(t=+.5))

* style and color for y-axis title/labels
local y1_title_style it
local y1_title_color gs8
local y1_label_style it
local y1_label_color gs8
local y1_label_size 9pt

local text_style sf
local text_color black

* y-axis labels placed above grid rule (build manually)
local y1_x `xlo'
local y1_place ne
local y1_just left

local y 1
local y`y'_la
foreach tick of numlist  `y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi' {
  if `tick' == `y`y'_tick_hi' {
    * white textbox to go behind label
    local y`y'_la `y`y'_la' text(`=`tick' + `y`y'_step'/100' `y`y'_x' " ", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') width(5) height(3) bcolor(white%70))
    *label
    local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`=string(`tick', "`y`y'_fmt'")'`y`y'_unit'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
  }
  else {
    local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`=string(`tick', "`y`y'_fmt'")'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
  }
}

* y-axis titles, horizontal orientation (build manually)
local y 1
local y1title text(`=`y`y'_tick_lo' + `y`y'_step'*(`y_ticks' - 1 + `y_margin_hi')' `xlo' "{`y`y'_title_style':Entrants per million residents}"
local y1title `y`y'title', yaxis(`y') place(`y`y'_place') just(`y`y'_just') box lwidth(none) color(`y`y'_title_color') bcolor(white) margin(b=+.5))
local y1title `y`y'title' ytitle(" ", axis(`y') margin(l=-8))

* y-axis settings
local y_axis1 `y1_label' `y1_la' `y1_scale' `y1title'

* Plot lines
local y1 incidence
local xvar `x_bin_mean'
local line_y1 (connected `y1' `xvar', pstyle(p2))
local ci_area (rarea `y1'_ci_lo `y1'_ci_hi `xvar', sort pstyle(p2area))

* Label line_y1
local text_1 text(362 8.5 "{`text_style':Monthly DI entry}"
local text_1 `text_1', yaxis(1) place(nw) just(right) box lwidth(none) width(30) color(`text_color') bcolor(white%90))
local xref = `xvar'[18]
qui sum `y1' if `xvar' == `xref'
local pci_1 (pci `=r(mean)' `=`xref'' `=r(mean) + 5' `=`xref' - 0.5', pstyle(p5))

local legend legend(off)
local region graphregion(margin(t=5 b=-1))

twoway `ci_area' `line_y1' `pci_1', `x_axis' `y_axis1' `legend' `region' `text_1'

* Export graph
cap mkdir "$SSDIMed/results/figures"
graph export "$SSDIMed/results/figures/unemp_ventiles_incidence-`sample'.pdf", as(pdf) fontface("Source Sans Pro") fontdir("$SSDIMed/data/raw/fonts") replace

}


* ------------------------------
* Panel (b): DI entrant cohort average spending (tot_pmt)
if 1 {

* How many major ticks?
* How much should range extend beyond ticks (as fraction of tick step size)
local y_ticks 7
local y_margin_lo 0.4
local y_margin_hi 0.4

* y1 tick lo, step size, format
local y1_tick_lo 13.0
local y1_step 0.1
local y1_fmt "%9.1f"
local y1_sign "$"
local y1_unit " thousand"

* Construct y labels and scales from parameters above
forvalues y = 1/1 {
  local y`y'_tick_hi  = `y`y'_tick_lo' + `y`y'_step' * (`y_ticks' - 1)
  local y`y'_range_lo = `y`y'_tick_lo' - `y`y'_step' * `y_margin_lo'
  local y`y'_range_hi = `y`y'_tick_hi' + `y`y'_step' * `y_margin_hi'
  di "y`y'_tick_hi:  `y`y'_tick_hi'"
  di "y`y'_range_lo: `y`y'_range_lo'"
  di "y`y'_range_hi: `y`y'_range_hi'"
  
  local y`y'_label ylabel(`y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi', axis(`y') format("`y`y'_fmt'") noticks nolabels labgap(1pt)) 
  local y`y'_scale yscale(`yalt' range(`y`y'_range_lo' `y`y'_range_hi') axis(`y') noline) 
  di `"y`y'_label:  `y`y'_label'"'
  di `"y`y'_scale:  `y`y'_scale'"'
}

* x-axis settings
local xlo 1.3
local xhi 13.8
local xscale xscale(range(`xlo' `xhi'))
local x_axis xlabel(3(1)13) `xscale' xtitle("County unemployment at application", margin(t=+.5))

* style and color for y-axis title/labels
local y1_title_style it
local y1_title_color gs8
local y1_label_style it
local y1_label_color gs8
local y1_label_size 9pt

local text_style sf
local text_color black

* y-axis labels placed above grid rule (build manually)
local y1_x `xlo'
local y1_place ne
local y1_just left
local y1_width 25

local y 1
local y`y'_la
foreach tick of numlist  `y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi' {
  if `tick' == `y`y'_tick_hi' {
    * white textbox to go behind label
    local y`y'_la `y`y'_la' text(`=`tick' + `y`y'_step'/100' `y`y'_x' " ", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') width(`y`y'_width') height(4) bcolor(white%70))
    *label
    local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`y`y'_sign'`=string(`tick', "`y`y'_fmt'")'`y`y'_unit'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
  }
  else {
    local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`y`y'_sign'`=string(`tick', "`y`y'_fmt'")'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
  }
}

* y-axis titles, horizontal orientation (build manually)
local y 1
local y1title text(`=`y`y'_tick_lo' + `y`y'_step'*(`y_ticks' - 1 + `y_margin_hi')' `xlo' "{`y`y'_title_style':Medical spending}"
local y1title `y`y'title', yaxis(`y') place(`y`y'_place') just(`y`y'_just') box lwidth(none) color(`y`y'_title_color') bcolor(white) margin(b=+.5))
local y1title `y`y'title' ytitle(" ", axis(`y') margin(l=-8))

* y-axis settings
local y_axis1 `y1_label' `y1_la' `y1_scale' `y1title'

* Plot lines
local y1 tot_pmt
local xvar `x_bin_mean'
local line_y1 (connected `y1' `xvar', pstyle(p3))
local ci_area (rarea `y1'_ci_lo `y1'_ci_hi `xvar', sort pstyle(p3area))

* Label line_y1
local text_1 text(13.51 5.51 "{`text_style':Annual medical spending}"
local text_1 `text_1', yaxis(1) place(ne) just(left) box lwidth(none) width(45) color(`text_color') bcolor(white%90))
local xref = `xvar'[8]
qui sum `y1' if `xvar' == `xref'
local pci_1 (pci `=r(mean)' `=`xref'' `=r(mean) + 0.07' `=`xref' + 0.45', pstyle(p5))

local legend legend(off)
local region graphregion(margin(t=5 b=-1))

twoway `ci_area' `line_y1' `pci_1', `x_axis' `y_axis1' `legend' `region' `text_1'

* Export graph
cap mkdir "$SSDIMed/results/figures"
graph export "$SSDIMed/results/figures/unemp_ventiles_tot_pmt-`sample'.pdf", as(pdf) fontface("Source Sans Pro") fontdir("$SSDIMed/data/raw/fonts") replace

}


* ------------------------------
* Panel (c): DI entrant cohort Medicare mortality (died)
if 1 {

* How many major ticks?
* How much should range extend beyond ticks (as fraction of tick step size)
local y_ticks 6
local y_margin_lo 0.33
local y_margin_hi 0.33

* y1 tick lo, step size, format
local y1_tick_lo 276
local y1_step 2
local y1_fmt "%9.2g"
local y1_unit ""

* Construct y labels and scales from parameters above
forvalues y = 1/1 {
  local y`y'_tick_hi  = `y`y'_tick_lo' + `y`y'_step' * (`y_ticks' - 1)
  local y`y'_range_lo = `y`y'_tick_lo' - `y`y'_step' * `y_margin_lo'
  local y`y'_range_hi = `y`y'_tick_hi' + `y`y'_step' * `y_margin_hi'
  di "y`y'_tick_hi:  `y`y'_tick_hi'"
  di "y`y'_range_lo: `y`y'_range_lo'"
  di "y`y'_range_hi: `y`y'_range_hi'"
  
  local y`y'_label ylabel(`y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi', axis(`y') format("`y`y'_fmt'") noticks nolabels labgap(1pt)) 
  local y`y'_scale yscale(`yalt' range(`y`y'_range_lo' `y`y'_range_hi') axis(`y') noline) 
  di `"y`y'_label:  `y`y'_label'"'
  di `"y`y'_scale:  `y`y'_scale'"'
}

* x-axis settings
local xlo 1.3
local xhi 13.8
local xscale xscale(range(`xlo' `xhi'))
local x_axis xlabel(3(1)13) `xscale' xtitle("County unemployment at application", margin(t=+.5))

* style and color for y-axis title/labels
local y1_title_style it
local y1_title_color gs8
local y1_label_style it
local y1_label_color gs8
local y1_label_size 9pt

local text_style sf
local text_color black

* y-axis labels placed above grid rule (build manually)
local y1_x `xlo'
local y1_place ne
local y1_just left

local y 1
local y`y'_la
foreach tick of numlist  `y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi' {
  if `tick' == `y`y'_tick_hi' {
    * white textbox to go behind label
    local y`y'_la `y`y'_la' text(`=`tick' + `y`y'_step'/100' `y`y'_x' " ", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') width(5) height(3) bcolor(white%70))
    *label
    local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`=string(`tick', "`y`y'_fmt'")'`y`y'_unit'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
  }
  else {
    local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`=string(`tick', "`y`y'_fmt'")'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
  }
}

* y-axis titles, horizontal orientation (build manually)
local y 1
local y1title text(`=`y`y'_tick_lo' + `y`y'_step'*(`y_ticks' - 1 + `y_margin_hi')' `xlo' "{`y`y'_title_style':Deaths per 10,000}"
local y1title `y`y'title', yaxis(`y') place(`y`y'_place') just(`y`y'_just') box lwidth(none) color(`y`y'_title_color') bcolor(white) margin(b=+.5))
local y1title `y`y'title' ytitle(" ", axis(`y') margin(l=-8))

* y-axis settings
local y_axis1 `y1_label' `y1_la' `y1_scale' `y1title'

* Plot lines
local y1 died_adj
local xvar `x_bin_mean'
local line_y1 (connected `y1' `xvar', pstyle(p4))
local ci_area (rarea `y1'_ci_lo `y1'_ci_hi `xvar', sort pstyle(p4area))

* Label line_y1
local text_1 text(283.0 5.25 "{`text_style':Annual mortality}"
local text_1 `text_1', yaxis(1) place(ne) just(left) box lwidth(none) width(35) color(`text_color') bcolor(white%90))
local xref = `xvar'[5]
qui sum `y1' if `xvar' == `xref'
local pci_1 (pci `=r(mean)' `=`xref'' `=r(mean) + 1.1' `=`xref' + 0.8', pstyle(p5))

local legend legend(off)
local region graphregion(margin(t=5 b=-1))

twoway `ci_area' `line_y1' `pci_1', `x_axis' `y_axis1' `legend' `region' `text_1'

* Export graph
cap mkdir "$SSDIMed/results/figures"
graph export "$SSDIMed/results/figures/unemp_ventiles_died-`sample'.pdf", as(pdf) fontface("Source Sans Pro") fontdir("$SSDIMed/data/raw/fonts") replace

}









** EOF
